import numpy as np
import newton as new
import math

print("PARTIE ELECTROSTATIQUE")

#---------Question 1-------#

#E(x1,...,xN) = somme_de_i=1_a_N(log|xi + 1| + log|xi - 1| + 1/2 * somme_de_j=1_a_N_j!=i(log|xi - xj|))
#Ei = dE/dxi = 1/(xi + 1) + 1/(xi - 1) + somme_de_j=1_a_N_j!=i(1/(xi - xj))
#Ci-dessous, la matrice E est la matrice jacobienne du vecteur Ei(1<i<N)
#pour i!=k, (pas les termes diagonaux), dEi/dxk = E[i,k] = 1/(xi - xk)**2
#pour i=k, (termes diagonaux), dEi/dxk = -(1/(xi + 1)**2 + 1/(xi - 1)**2 + somme_de_j=1_a_N_j!=i(1/(xi - xj)**2))

N=3#a modifier librement selon le nombre de variables

def somme_termes_diagonaux(x,i):#somme_de_j=1_a_N_j!=i(1/(xi - xj)**2)
    s = 0
    for j in np.arange(0,N):
        if i != j:
            s = s + 1.0/(x[i]-x[j])**2
    return s
    
def electrostatic_energy_jacobian(x):#retourne la matrice jacobienne de nablaE(x), x vecteur a N coordonnees
    t = np.arange(0,N)
    E = np.zeros([N,N])
    for i in t:
        for j in t:
            if i !=j:
                E[i,j] = 1/(x[i]-x[j])**2
    for i in t:
        E[i,i] = -(1.0/(x[i]+1)**2 + 1.0/(x[i]-1)**2 + somme_termes_diagonaux(x,i))
    return E

print electrostatic_energy_jacobian(np.array([0.2,0.4,0.5]))

#------------ Question 2 --------------#

def somme_termes_derivee(x,i):#somme_de_j=1_a_N_j!=i(1/(xi - xj))
    s = 0
    for j in np.arange(0,N):
        if i != j:
            s = s + 1.0/(x[i]-x[j])
    return s

def somme_termes_log(x,i):#somme_de_j=1_a_N_j!=i(log|xi - xj|)
    s = 0
    for j in np.arange(0,N):
        if i != j:
            s = s + math.log(math.fabs(x[i]-x[j]))
    return s

def electrostatic_energy_nabla(x):#retourne les dE/dxi sous forme d'un tableau
    t = np.arange(0,N)
    E = np.zeros(N)
    for i in np.arange(0,N):
        E[i] = 1.0/(x[i]+1) + 1.0/(x[i]-1) + somme_termes_derivee(x,i)
    return E

def electrostatic_energy_total(x):#retourne E(x1,...,xN), l'energie totale du systeme
    t = np.arange(0,N)
    s = 0
    for i in t:
        s = s + math.log(math.fabs(x[i]+1)) + math.log(math.fabs(x[i]-1)) + 1.0/2 * somme_termes_log(x,i)
    return s

print electrostatic_energy_nabla(np.array([0.2,0.4,0.5]))

print ("Resolution de l'equation :")

print new.Newton_Raphson_N(electrostatic_energy_nabla,electrostatic_energy_jacobian,np.array([0.2,0.4,0.5]),0.1,50)

#on constate (en regardant la page wiki par exemple) que les solutions correspondent bien aux racines des polynomes de Legendre

#------Question 3-----#

print("Maximum ou minimum :")

ref = new.Newton_Raphson_N(electrostatic_energy_nabla,electrostatic_energy_jacobian,np.array([0.2,0.4,0.5]),0.1,50)
test = np.array([0.2,0.4,0.5])
test2 = np.array([0.1,-0.4,-0.5])


print electrostatic_energy_total(ref)
print electrostatic_energy_total(test)
print electrostatic_energy_total(test2))

M = electrostatic_energy_jacobian(ref)
print M
print np.linalg.eigvals(M)#calcule les valeurs propres de M
